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We extend White’s minimally entangled typically thermal states approach (METTS) to allow 
Abelian and non-Ablian symmetries to be exploited when computing finite-temperature response 
functions in one-dimensional (ID) quantum systems. Our approach, called SYMETTS, starts from 
a METTS sample of states that are not symmetry eigenstates, and generates from each a sym¬ 
metry eigenstate. These symmetry states are then used to calculate dynamic response functions. 
SYMETTS is ideally suited to determine the low-temperature spectra of ID quantum systems with 
high resolution. We employ this method to study a generalized diamond chain model for the natural 
mineral azurite 0113(003)2(011)2, which features a plateau at | in the magnetization curve at low 
temperatures. Our calculations provide new insight into the effects of temperature on magnetiza¬ 
tion and excitation spectra in the plateau phase, which can be fully understood in terms of the 
microscopic model. 

PACS numbers: 71.27.+a, 75.10.Pq 


I. INTRODUCTION 

The simulation of dynamical quantities in one¬ 
dimensional (ID) quantum many-body systems still 
poses a major challenge for theoretical condensed mat¬ 
ter physics, particularly at finite temperature. From an 
experimentalist’s perspective, there is high demand for 
such calculations for a variety of reasons: (i) Experimen¬ 
tal measurements hardly allow to study solely ground- 
state physics as thermal fluctuations cannot be elimi¬ 
nated altogether. Thus for a direct comparison with ex¬ 
perimental data, it is essential to include temperature in 
the theoretical modeling, (ii) Technical advances have 
nowadays drastically enhanced the precision of neutron 
scattering and electron resonance spectroscopy, which 
for example allows the measurement of dynamic observ¬ 
ables such as momentum-resolved excitation spectra in 
effective ID materials with very high resolution. 1-6 (iii) 
Thermal fluctuations can cause new phenomena, which 
are not captured by the ground-state physics of the sys¬ 
tem. Two examples are the sudden emergence of a sin¬ 
gle spinon dispersion (“Villain mode”) in XXZ-like spin- 
chain materials, 7-9 or the existence of quantum critical 
phases in various strongly correlated materials. 10 

Which numerical tools can be employed to simulate 
such dynamic observables in a ID quantum system? 
At zero temperature, the density matrix renormaliza¬ 
tion group (DMRG) is the most successful exact nu¬ 
merical method for describing quantum many-body sys¬ 
tems regarding their static and dynamic ground-state 
properties. 11,12 DMRG based algorithms have also been 
successfully extended to treat systems at finite tempera¬ 
ture, yet the computational efficiency of such approaches 
is still limited. Exact diagonalization (ED) or quantum 
Monte Carlo (QMC) 13-15 can only be considered as com¬ 
plementary approaches rather than proper alternatives 
to DMRG, since the applicability of ED is restricted by 
small system sizes, and that of QMC by the need for 


performing an ill-defined analytic continuation, and of¬ 
ten also by the occurrence of a sign problem. Thus, the 
simulation of experimentally relevant quantities such as 
dynamic response functions represents a highly demand¬ 
ing and difficult task for finite-temperature numerics. 

Whereas early DMRG approaches for computing finite- 
temperature response functions for a ID quantum sys¬ 
tem have been based on the transfer matrix renormal¬ 
ization group (TMRG), 16-19 today the most popular 
method builds on the purification of the density matrix 
in the matrix-product-state (MPS) formalism. 20 The re¬ 
sponse functions can then be calculated with high pre¬ 
cision by using tDMRG in the real-time realm, and a 
subsequent Fourier transform also allows the computa¬ 
tion of spectral functions. 21-2 ' Purification can also be 
combined with a Chebyshev expansion technique to de¬ 
termine finite-temperature spectral functions directly in 
frequency space. 28,29 Although these methods have been 
successfully applied to a number of experimental setups, 
the accessible time scale (or maximal Chebyshev expan¬ 
sion order) and hence the spectral resolution is limited, as 
the propagation of excitations during the dynamic evolu¬ 
tion yields a linear growth of entanglement, leading to an 
exponential increase in the required numerical resources. 
In addition, the encoding of mixed states inevitably re¬ 
quires doubling the size of the Hilbert space and intro¬ 
duces additional entanglement between the physical state 
and its environment, which limits the efficiency of purifi¬ 
cation simulations towards low temperatures. 

An alternative way to compute finite-temperature 
quantities was recently presented by White in Ref. [30]. 
Instead of purifying the density matrix, an ensemble of 
pure states is introduced that are constructed to resem¬ 
ble the typical state of a quantum system at finite tem¬ 
perature. It has been shown that these so-called mini¬ 
mally entangled typical thermal states (METTS) excel¬ 
lently represent the thermal properties of the system of 
interest. At the same time, they can efficiently be rep- 
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resented in the MPS formalism as their entanglement is 
very low. 30,31 The METTS approach was originally only 
used to compute static quantities of spin-chains 30,31 and 
fermions. 32 In the meantime, it has also been applied 
to simulate finite-temperature quenches 33 and response 
functions. 34 

The numerical effort for constructing a single METTS 
is comparable to ground-state DMRG, since METTS 
avoids the explicit computation of the density matrix. 
Since METTS calculations are also easily parallelized, 
it has originally been considered to be a more efficient 
finite-temperature formulation than purification. More 
recently, Ref. [34] showed in a detailed study that this 
claim cannot generally be supported, because the addi¬ 
tional statistical error source introduced by the sampling 
increases computational costs, especially at high temper¬ 
atures. Nevertheless, METTS still offers much potential 
towards the simulations of low-temperature properties of 
complex models, as long as one does not insist on reduc¬ 
ing the statistical error to be as small as the truncation 
error. 

To bring out the full potential of METTS for the calcu¬ 
lation of dynamic quantities, this work addresses a severe 
constraint of the current formulation of the algorithm: 
the ensemble states cannot be chosen such that they re¬ 
spect inherent symmetries of the system and at the same 
time minimize autocorrelation effects. This drastically 
increases the numerical resources necessary for comput¬ 
ing the real-time evolution of the ensemble states, as the 
MPS have not been decomposed into symmetry blocks by 
means of the symmetry-induced selection rules. 35-38 To 
remedy this problem, we introduce an intuitive and easily 
implementable extension of White’s approach: starting 
from a METTS sample of states that are not symme¬ 
try eigenstates, we generate a sample of symmetry eigen¬ 
states, called SYMETTS. These states allow both simple 
Abelian and more complex non-Abelian symmetries to 
be exploited in the computation of dynamic quantities. 

As an experimentally relevant application of 
SYMETTS, we study temperature effects on the | 
magnetization plateau of a generalized diamond chain 
model, which has been derived as a microscopic model 
for the natural mineral azurite Cu 3 (C 03 ) 2 ( 0 H ) 2 . 39,40 
This material has attracted much attention due to the 
discovery of a plateau at | in the magnetization curve 
at low temperatures. 4,5,39-49 Via real-time evolution of 
SYMETTS ensembles, it is possible to obtain highly 
resolved excitation spectra in the ^ plateau phase for 
various temperatures. We observe a crossing of monomer 
and dimer branches with increasing magnetic field, which 
intuitively explains the effects of finite temperature on 
the magnetization in the plateau phase. 

The paper is organized as follows. In Sec. II we briefly 
review the original METTS algorithm and the neces¬ 
sity of choosing a symmetry-breaking collapse routine to 
generate the ensemble. Sec. Ill introduces a METTS 
formulation based on symmetry eigenstates for models 
with both Abelian and non-Abelian symmetries. Sec. IV 


summarizes benchmark calculations for static and dy¬ 
namic observables for the spin-^ XXZ chain. In Sec. V 
SYMETTS is employed to study an experimentally rele¬ 
vant microscopic model for the natural mineral azurite. 
A technical discussion on the combination of SYMETTS 
with a Chebyshev expansion to directly calculate dy¬ 
namic correlators in frequency space is relegated to ap¬ 
pendix A. The computational efficiency of SYMETTS in 
the context azurite is assessed in appendix B. 

II. MINIMALLY ENTANGLED TYPICAL 
THERMAL STATES 

A. METTS calculations thermal quantities 

First of all, we review the construction of a METTS 
sample to approximate a thermal expectation value (A) p 
for a general chain model with N sites. To this end, 
the trace of a thermal expectation value (A)p = Tr [ppA\ 
is expanded in terms of an orthonormal basis {| cr)} 
of classical product states (CPS) of the form |<x) = 
|cr 1 ) \a 2 )...\a N ). Each such state has an entanglement 
entropy of exactly zero. Thus these states represent a 
natural choice for a basis at infinite temperature, where 
the system should behave classically. In addition, their 
entanglement growth under imaginary-time evolution re¬ 
mains comparatively low, hence the designation “mini¬ 
mally entangled” states. The expectation value of A can 
be written as 

(A) P = ^Y,^\ e ~^ /2Ae ~ 06/ » 

= ^~'E P ^\A\K), (!) 

with the partition function Zp = Tr[e _/3i? ] = Y2o-P<r- 
The normalized states \<f>a) represent a set of METTS 
with corresponding probabilities P rTl defined as 

l<M = ^=e -,3ff/2 |cr), P rT = {a\e- pH \cr). (2) 

By sampling the METTS \(j) CT } according to the proba¬ 
bility distribution P^/Zp, the calculation of a thermal 
expectation value can be reformulated into taking the 
plain average of (cj) rT \A\(t> CT ). 

To obtain a METTS sample {\4>a)} with the correct 
probability distribution, a Markov chain of CPS | cr) is 
generated. This is done in a way that obeys detailed 
balance, which guarantees reproducing the probability 
distribution P a /Zp. The sampling algorithm can be set 
up sequentially. To this end, one starts from an arbi¬ 
trary CPS |cr) and conducts what is called a thermal step: 

(i) A single METTS |^> CT ) is generated by evolving the 
CPS in imaginary time and normalizing it. 

(ii) A measurement of all local degrees of freedom is 
performed by projecting (or collapsing) \<j>a) into a 



3 


k) ( ii) k'> k") 

(i)e -^ /2 P^e-¥' 2 A 6 * 

I / I / I / 

I |0i t) |0<x') |0ct") I 

--J-- 

(0<r|^|0<r) *- {I0<x)} -► <0« T |S(t)C'|0, T ) 

Figure 1. Schematic illustration of the METTS algorithm. 
For details on how to explicitly evaluate response functions of 
the type {B(t)C)p see Sec. IV B. 

new CPS | er') with probability p CT ' CT = \{cr l \(j ) rT )| 2 . 

The transition probabilities obey detailed balance 

Pcr'a-Pcr = Pcra'Pcr' by Construction. 

The thermal step is then repeated with the newly gen¬ 
erated CPS to generate a METTS (see Fig. 1 for illustra¬ 
tion). By construction, the correct distribution is recov¬ 
ered as a fixed point of this procedure. To eliminate any 
artificial bias caused by the choice of the initial random 
CPS, the first few thermal steps are neglected in the cal¬ 
culation of any static observable (<^ er |.A|0 (T ) or dynamic 
response function (B(t)C)p. 

By making good choices for the local measurements 
(see Sec. IIB), the sample size M can be chosen surpris¬ 
ingly small to obtain accurate results [M ~ 0(1O 2 -1O 3 )]. 

At this point, a comment on accuracy is in order. It 
was recently pointed out by Ref. [34] that the purifica¬ 
tion approach drastically outperforms METTS, because 
for fixed total computation time it reaches more accurate 
results, where the accuracy was judged by comparing to 
quasi-exact calculations. However, this should not be a 
surprise, since the METTS sampling introduces an addi¬ 
tional statistical error source, which generally scales as 
f3~ 1 /\/M — 1. Obviously, this prevents a perfect conver¬ 
gence of METTS results towards exact data and limits 
the efficiency at very high temperatures in comparison to 
a non-statistical method. Nevertheless, we believe that 
METTS offers much potential towards the simulations of 
low-tenrperature properties of complex models, as long 
as one does not insist on pushing the statistical error 
towards the order of the truncation error. 


B. Ergodicity and efficient sampling 

Generating a new CPS |er) by collapsing a METTS 
represents the most crucial step of the sampling algo¬ 
rithms, as a bad choice of measurement basis leads to a 
drastically increased autocorrelation time. 31 

Let us assume that the local Hilbert space of each site 
j in our chain model is represented by an orthonormal 
basis | <jj) of size d, Uj € {1, 2,..., d}. The projective mea¬ 
surement 100 -) —> | <j') can be efficiently carried out site 
by site by making use of the well-defined orthogonality 


relations for a MPS, typically starting at one end of the 
chain (in our case site 1). To this end, the d transition 
probabilitiesp(<Ti) = (0 o .|P(<ti)| 0 o .) are calculated by in¬ 
troducing the projectors P(a i) = |cn)ki|- Then one of 
the d states is chosen with probability p{a\) by rolling 
a dice. The state is collapsed by the application of the 
projector P(<ti), and the orthonormal center of the MPS 
is shifted to the next site, where the collapse process is 
repeated. 

In principle, the orthonormal basis \aj) on each site j 
can be chosen arbitrarily. Nevertheless, there are good 
and bad choices with respect to the sampling efficiency. 
We illustrate this for the example of the spin-1 XXZ 
Heisenberg chain, 

N N 

h = JY. [sprt.Ws“«+ aS Pj«]+'‘X s J • P) 

3 j 

for the isotropic case J = 1, A = 1 and h = 0. This model 
features a non-Abelian SU(2) sp i n symmetry, which can be 
reduced to an Abelian U(l) symmetry, e.g., by consider¬ 
ing the total magnetization S^ ot as a good quantum num¬ 
ber. At first sight, the eigenstates of the spin operator 
Sj resemble a natural choice for the orthonormal basis 
set tjj ), since this choice allows the encoding of the pro¬ 
jectors in the form of diagonal operators. Moreover, all 
resulting CPS are eigenstates of S z . There¬ 

fore, it is possible to directly implement the Abelian U(l) 
symmetry in the MPS representation resulting in a mas¬ 
sive reduction of computational effort. 

However, a collapse routine based on measurements 
along the 2 axis only (“z collapse”) leads to a serious 
problem with ergodicity, as already extensively discussed 
in Ref. [31]. Subsequently generated CPS are strongly 
correlated and thus the autocorrelation times are very 
long, so that the bias arising from the initial random 
CPS cannot be removed in a few thermal steps. Addi¬ 
tionally, CPS generated from subsequent thermal steps 
always have the same total magnetization S^ ot since the 
z collapse conserves this quantity. It is therefore im¬ 
possible to cover different S z ot sectors with the sampling 
algorithm described above, as one is always stuck in the 
symmetry sector of the initially chosen CPS. 

This issue can be resolved by randomly choosing a dif¬ 
ferent local basis for each site of the chain (“random col¬ 
lapse”). Alternatively, it has been shown that alternat¬ 
ing in subsequent thermal steps between two basis sets 
that are maximally mixed relative to each other, e.g., the 
eigenstates of S z and S x , also restores ergodicity and cov¬ 
ers multiple symmetry sectors of the sample (“maximally 
mixed collapse”). 

We illustrate the failure of the 2 collapse routine by 
starting from an ensemble of randomly generated CPS 
and then conducting 10 thermal steps with each state 
of the ensemble for three different values of /3. After 
each step, we measure the ensemble average of the total 
energy, which is displayed in Fig. 2. For comparison, 
we also calculate the same quantity using the random 
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dynamic properties of complex models, such as ID sys¬ 
tems and 2D lattice models with experimental relevance. 
Since the current METTS setup does not allow Abelian 
or non-Abelian symmetries to be exploited, it is not suit¬ 
able for accessing dynamic observables of such complex 
systems. 

In the following, we show how to resolve this funda¬ 
mental issue by a simple extension of the sampling al¬ 
gorithm that will enable us to systematically build a 
METTS ensemble based on symmetry eigenstates. 


III. SYMMETRIC METTS 


Figure 2. Energy of the spin-Heisenberg chain with N = 
100 spins for /3 = 4,8,20. Details on the setup of the 
imaginary-time evolution can be found in Sec. IV. Starting 
from an ensemble of 100 randomly generated CPS, we con¬ 
duct ten thermal steps with each state and measure the en¬ 
semble average of the total energy in each step. The different 
basis choices for the CPS collapse become apparent in the 
autocorrelation times. Whereas measuring along the z axis 
only (squares) leads to strong autocorrelations that prevent 
the energy to converge towards the exact value (black lines), 
randomly chosen measurement bases (circles) result in short 
autocorrelation times of a few thermal steps. 31 


collapse routine. When choosing a random basis for each 
CPS collapse (circles), the total energy of the ensemble 
is already well converged towards the exact value after 
a few thermal steps, because autocorrelations between 
subsequent CPS are practically absent. 

In contrast, when measuring along the z axis only 
(squares), the total energy is nowhere near its exact 
value, even after 10 thermal steps. We have discussed the 
causes of this behavior above: first of all, one can iden¬ 
tify strong correlations between subsequent CPS during 
the application of the z collapse resulting in an increase 
of autocorrelation time. In addition, each CPS remains 
in its initial symmetry sector. If the different symmetry 
sectors are not distributed according to the correct prob¬ 
ability distribution at a specific value /3 (which is very 
unlikely starting from a random set), the ensemble can¬ 
not capture the correct behavior of the system, as the 
sample is biased towards specific sectors. This explains 
why the average energy is not only converging slowly to¬ 
wards the exact value, but rather seems to saturate at a 
significantly higher value. Thus, a symmetry conserving 
collapse routine that is based on measurements in a fixed 
local basis is clearly impracticable. 

If we want to retain the ergodicity of the METTS 
sample, we are left to choose between the random or 
the maximally mixed collapse routine. This comes at 
a price, as the ensemble states cannot be chosen such 
that they conserve inherent symmetries of the system, 
because both collapse routines clearly require symmetry¬ 
breaking measurements. However, the efficient treatment 
of symmetries is often essential for calculating especially 


A. Symmetries 

The matrix-product-state framework allows for a 
straightforward incorporation of symmetries of the model 
Hamiltonian. 35-38 Generally speaking, the symmetry- 
induced selection rules cause a large number of matrix 
elements to be exactly zero, thus bringing the Hamilto¬ 
nian into a block-diagonal structure and subdividing ten¬ 
sors into well-defined symmetry sectors. Keeping only 
the non-zero elements, we can achieve tremendous im¬ 
provement in speed and accuracy in numerical simula¬ 
tions by the inclusion of symmetries. In the context of 
non-Abclian symmetries, the non-zero data blocks are 
not independent of each other and can be further com¬ 
pressed using the Clebsch-Gordan algebra for multiplet 
spaces. Here we refrain from discussing this topic at 
length and refer to Ref. [36] for a detailed review on the 
treatment of symmetries in tensor network applications. 

Following the notation of Ref. [36] , we label the state 
space in terms of the symmetry eigenbasis | qn',q z ), where 
the quantum labels q denotes the irreducible representa¬ 
tion of the symmetry group S of the Hamiltonian H. Ev¬ 
ery symmetry generator S a satisfies [H, 5 a ] = 0 . Hence, 
all states in a given Hilbert space corresponding to a 
certain g-label are combined into a symmetry block q. 
The label n identifies a particular multiplet within the 
specific symmetry block q. The internal multiplet la¬ 
bel q z resolves the internal structure of the correspond¬ 
ing multiplet. In the context of Abelian symmetries the 
Clebsch-Gordan structure becomes trivial, hence the q z 
labels take the role of q labels. Note that this notation 
can be easily generalized to the treatment of multiple 
symmetries. 36 

To further clarify the notation, we consider the exam¬ 
ple of the isotropic Heisenberg chain in Eq. (3), which 
features a SU(2) sp ; n symmetry, S = SU(2) sp i n . We make 
the usual choice of basis in which the 2 component of 
the spin operator, S z is diagonal and label a general spin 
multiplet by \q,q z ) = |S, S z ). The spin multiplet label 
can take the values q = 0, 4,1, |, ..., while the internal 
multiplet label, corresponding to the z component of the 
spin, is restricted to q z £ {— q, —q + 1 ,..., +q}. 

Now consider a typical MPS scenario, where the wave 
function \ip) in the local picture of site j can be repre- 
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IV>) = E ^\L)\aj)\R). (4) 

LcTjR 

In the presence of symmetries, the physical state space 
at site j as well as the left and right orthonormal basis 
states can be written as | L) = | ql;q z ), |oy) = \q'in;q' z ), 
|i?) = \q"n\q"). Hence, symmetry labels can be intro¬ 
duced naturally in the MPS representation. In particu¬ 
lar, every leg or bond in the usual diagrammatic depiction 
of a MPS can be assigned a multiplet label, here 9 , q' and 
<l", e.g., 



B. METTS with symmetry eigenstates 


In order to work with a symmetry conserving METTS 
ensemble, we reformulate Eqs. (1) and (2) in terms of 
symmetry eigenstates before introducing an efficient sam¬ 
pling routine (see Sec. IIIC). In place of the CPS we in¬ 
troduce a set of symmetry product states (SPS) | q), that 
can be considered as symmetrized counterparts of the 
CPS. A SPS is a MPS with (multiplet) bond dimension 
one, where each bond represents a single, unique symme¬ 
try block < 7 j-, and that block contains just a single multi¬ 
plet (rij = 1). Thus the SPS can be fully characterized 
by a set of N quantum labels q = { q± , < 72 , 93 ; •••, < 7 w}, one 
label qj per site/bond j labeling the corresponding sym¬ 
metry sector. The overall symmetry sector of each SPS 
| q) is fully determined by the 9 -label of the last bond, 
9.v- 

The simplest example of a SPS for the SU(2) symmet¬ 
ric Heisenberg chain is to combine pairs of neighboring 
spins to singlets 

Iq)sU( 2) = (I tl) | 4-2 ) - | 1*2)1 4'l))(| t3)l -U) - 11 3 ) I I 4 )) 

(I tiv-i)l In) - I tjv— 1 )| In))- (6) 

For spin-1 systems, the quantum labels q correspond to 
a sequence of qj = ^ , 0 for odd and even bonds, respec¬ 
tively, with a total spin 9 yr = S'tot = 0 , as illustrated in 
Fig. 3(a). Although the (multiplet) dimension on each 
bond remains one for a non-Abelian SPS, it is no longer 
a pure product state in the classical sense as the internal 
multiplet structures can introduce non-trivial entangle¬ 
ment between neighboring sites. 

The same formalism also applies for Abelian SPS. Re¬ 
ducing SU(2) to an Abelian U(l) symmetry in the Heisen¬ 
berg model, e.g., by choosing an anisotropy A / 1 in 
Eq. (3) or adding a finite magnetic field in the z direction, 
and choosing the total spin S* ot as conserved quantity, a 
typical SPS takes the form of 


(7) 


l9)su(2) 

l«)u(i) 


(a) 91 = x 12 = 0 ® = x 


In = 0 


T T T T-T 


<7i - 12 = -1 93 = 


Qn = 2 


(b) 


T T T T-T 


Figure 3. Schematic illustration of the SPS (a) in Eq. (6) and 
(b) in Eq. (7). 


In this case, the quantum label qj represents the sum of 
all 5 z -contributions for sites i < j, i.e., q 3 = Sf as 
shown in Fig. 3(b). Hence, the total magnetization S^ ot 
of the SPS is given by the last label 9 jy. In the Abelian 
case, a SPS can be understood as a direct product of local 
symmetry eigenstates of each site, and hence is always 
represented by a MPS of bond dimension one, much like 
a classical product state. 

Analogously to the CPS basis set {| cr) }, a full set of 
SPS (|q)} represents a complete orthonormal basis tak¬ 
ing into account all possible symmetry sectors of the sys¬ 
tem. Thus, we can proceed as above and expand the 
trace of a thermal expectation value (A)p = Tr [/ 2 / 9 A] in 
terms of the symmetry product states | q): 


(A)p 


Zii q 


( 8 ) 


The normalized states \4> q ) now represent a set of sym¬ 
metric METTS (SYMETTS) with probabilities P q de¬ 
fined in analogy to Eq. (2), 

l<Aq) = —(9a) 
P q = (q\e-^\q) . (9b) 

The thermal expectation value (A)p is now estimated 
by sampling SYMETTS \(f> q ) according to the probability 
distribution P q /Zp. However, we still have to establish 
how to sample a set of SYMETTS {l^q)} according to 
the correct probability distribution Pq/Z^, in a way that 
ensures ergodicity. 


C. Algorithm for efficient sampling 

We illustrated in Sec. IIB for the spin-| Heisenberg 
chain that a collapse routine purely based on measure¬ 
ments along the z axis, conserving the U(l) sp i n sym¬ 
metry, fails to capture the correct thermal properties 
of the model. This is due to strong autocorrelation ef¬ 
fects and the fact that the symmetry sectors initially are 
distributed randomly instead of according to the correct 


l<?)u(l) = I 4-l)l 1-2)1 T 3) | -i-4) - - -1 Ihv)- 
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principle, it is possible to reduce the number of non- 
symmetric METTS by generating a larger number of 
SPS from each |^> CT ) by repeating the symmetrization 
step multiple times. However, this may introduce arti¬ 
ficial correlations between different SPS generated from 
the same non-symmetric METTS. For this reason, we 
present the SYMETTS algorithm using a formulation 
that maximizes ergodicity. This limits the applicabil¬ 
ity of SYMETTS in terms of calculating static observ¬ 
ables. In these cases, we would have to work harder 
than with regular METTS. However, when calculating 
dynamic quantities, generating the SYMETTS sample 
accounts for only a factor 0 ( 1 O~ 3 ) or less of the total 
computation time. Hence, our algorithm ensures that 
the full potential of SYMETTS towards dynamic appli¬ 
cations is guaranteed and no ergodicity problems arise. 


Figure 4. Schematic illustration of the SYMETTS sampling 
algorithm. For details on how to explicitly evaluate response 
functions of the type ( B(t)C)p see Sec. IV B. 

probability distribution P q /Zp. From this discussion, we 
can learn that the SYMETTS sample has to be generated 
from SPS that already capture the correct distribution of 
symmetry sectors . 50 

This can be achieved by starting from a typical ther¬ 
mal state, which already incorporates all the necessary 
thermal information. To this end, we extend the METTS 
sampling algorithm. After the conduction of a thermal 
step with a non-symmetric CPS |er) and METTS Ic/v) 
using a random or maximally mixed collapse [c.f (i) and 

(ii) in Sec. II], we employ an additional symmetrization 
step: 

(iii) Using a symmetry-conserving collapse routine (de¬ 
scribed below), we collapse | cj),?) to a SPS | q) with 
probability p qrT = |(q|<)> cr )| 2 . Each collapse gener¬ 
ates a SPS according to the correct probability dis¬ 
tribution P q /Zp (thus belonging to one of the rele¬ 
vant symmetry sectors at a given temperature), as 
long as the non-symmetric METTS has been sam¬ 
pled according to P^/Zp. 

(iv) The resulting SPS | q) can easily be converted 
into an MPS with explicit encoded symmetry 
sectors , 36 which is then evolved in imaginary time 
and normalized to generate the SYMETTS \4> q ). 

The combination of thermal and symmetrization step is 
then repeated with the newly generated CPS |er') to cre¬ 
ate a full SYMETTS sample {}, which represents the 
basis for calculating static or dynamic observables at fi¬ 
nite temperature (see Fig. 4 for an illustration). Thus we 
ensure that all computed SYMETTS are minimally auto- 
correlated, as each of them is generated from a different 
non-symmetric METTS. 

By maximizing the ergodicity of the sample, we face 
additional computational cost, as we have to generate a 
full non-symmetric METTS sample {|</> CT )} as well. In 


D. Collapse routine for non-Abelian symmetries 

In this section, we illustrate step (iii) of the SYMETTS 
sampling for the example of the isotropic spin-1 Heisen¬ 
berg chain (3). 

In context of U(l) sp i„, the collapse routine employed in 
step (iii) simply corresponds to the z collapse discussed 
in Sec. IIB, i.e., measuring along the z axis only. The 
resulting SPS take the form of direct products of local 
symmetry eigenstate | t) ? I 4-) and are automatically dis¬ 
tributed according to the correct probability P q /Zp. 

However, to exploit the full SU(2) spin symmetry of the 
model, the collapse routine has to be adapted in order to 
generate a SPS ensemble {|q)} of SU(2) eigenstates. For 
a single SPS, this is achieved by using a non-symmetric 
METTS \4>cr) and sequentially collapsing it into the dif¬ 
ferent eigensectors of the total spin operator S 2 . 

To this end, we gradually build S 2 starting from the 
left end of the chain. At the first site, the total spin 
is always 1 as we only consider a single spin, hence af¬ 
ter constructing S± no projection is required and the or¬ 
thonormal center of | </> CT ) can be shifted to the second 
site. Here, we generate the total spin operator of first 
and second site according to 

££,2 = s? + si + s x s 2 + s 2 Si, (io) 

with S 2 = (S'*) 2 + (Sj) 2 + (S~) 2 and the subscript “A, 2” 
indicating that we consider the total spin of the left part 
of the chain up to the second site. Diagonalizing this op¬ 
erator, we obtain the two spin sectors Sl, 2 = £ = 0,1 
corresponding to the singlet and triplet configuration, 
and the projectors P(Sl j2 )- We project the second bond 
of l^o-) (and also S 2 2 ) either into singlet or triplet con¬ 
figuration according to the transition probabilities 

p(Sl ti ) = (ii) 

and shift the orthonormal center of \4>a) to the next site. 
This procedure is repeated sequentially for every site j 
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of the system. Each time, we construct the spin operator 
for the left and the local part of the chain according to 

S 2 L j = Si d _, + S] + Slj-iSj + SjSLj-u ( 12 ) 

where S\ J _ 1 denotes the total spin squared of all sites to 
the left of (and excluding) site j. After cliagonalization, 
the transition probabilities are calculated and \4>a) is pro¬ 
jected at bond j into a single spin sector. Just as for the 
initially considered example, the operator S\ j always 
contains only two spin sectors, namely, Slj = i±|. 

Hence, diagonalization and projections can be carried out 
very efficiently. 

In the end, one obtains a SU(2) sp i n symmetric SPS | q) 
with probability p qrT = |(q|()> cr )| 2 . States of this type are 
the initial point for setting up the SU(2) symmetric MPS 
framework. 36 


IV. BENCHMARK RESULTS 

In this section, we present some benchmark results for 
our SYMETTS approach applied to both static and dy¬ 
namic observables of the XXZ Heisenberg model with 
N = 100 spins in the isotropic (A = 1, XXX model) and 
the free fermion limit (A = 0, XX model). As trunca¬ 
tion criterion, we choose to keep all singular values above 
sjj o1 > 10 -5 during the process of imaginary-time evolu¬ 
tion, which is carried out using standard tDMRG tools 
with a second order Trotter decomposition and a time- 
step t = 0.05. For the subsequent real-time evolution we 
adapt only the truncation criterion to = 10~ 4 . All 
quantities are expressed in terms of the coupling J = 1. 

A. Static observables: Thermal energy 

First, we discuss some static SYMETTS calculations 
for the total energy of the isotropic Heisenberg chain (3) 
with and without finite magnetic field. The data below 
conclusively show that the slightly modified METTS al¬ 
gorithm above is able to obtain results of similar accuracy 
as the non-symmetric METTS sampling at equal sample 
size M. Of course, this is to be expected since SYMETTS 
essentially generates the ensemble states analogously to 
the original algorithm. Nevertheless, this exercise helps 
to understand the importance of using sample states 
which are correctly distributed over the relevant sym¬ 
metry sectors. 

To illustrate the method in more detail, Fig. 5 shows 
U(1)-SYMETTS results resolving the different symme¬ 
try sectors entering into the calculation of the thermal 
energy for four different inverse temperatures. The up¬ 
per row displays the average energy of each subsample of 
states characterized by fixed S* ot . In the middle row, 
we zoom into a window of order of the temperature 
around the average energy of the sample. The resulting 
values for (E)p determined by SYMETTS [Eq. (8)] are 


(3=4 (3=8 (3=12 (3=20 



Figure 5. U(1)-SYMETTS sampling of the thermal energy 
for the isotropic Heisenberg chain with h = 1. Panels (a-d) in 
the upper row display the thermal energy of each symmetry 
sector entering the SYMETTS sample for /3 ranging between 
4 and 20; dashed lines indicate the overall ensemble average 
{E)p determined by SYMETTS. Moreover, a comparison to 
benchmark calculation based on METTS (dotted lines) and 
purification (solid lines) is provided. Panels (e-h) in the sec¬ 
ond row show /3(E — (E)p), in order to zoom into an energy 
window of order of the temperature around ( E)p. Panels (i-1) 
in the last row illustrates the subsample size M s ‘ ot of different 
symmetry sectors for a fixed total sampling size M = 500. 


benchmarked against METTS calculations [Eq. (1)] and 
quasi-exact purification data. The thermal average of all 
SYMETTS subsamples leads to highly accurate results 
for (E)p. 

For large f3, the lowest energy is obtained by the 
5* ot = 12 sector, which corresponds to the ground state 
symmetry sector of the system. Depending on the tem¬ 
perature, the energy of the neighboring sectors increases 
more or less steeply. At high temperatures, thermal 
fluctuations become clearly visible in the thermal en¬ 
ergies of the different symmetry sectors. Accordingly, 
the number of relevant symmetry sectors obtained from 
the SYMETTS sampling step (ii) is closely related to 
the temperature. Thermal fluctuations drive the sample 
states into more “excited” symmetry sectors at high tem¬ 
peratures: the maximum symmetry sector occurring in 
the /? = 4 simulation corresponds to S£ ot = 20, whereas 
we find a maximum of *S'^ ot = 15 for /3 = 20 as the system 
relaxes more towards the ground state. This behavior is 
also illustrated by the bottom row of Fig. 5, which shows 
the subsample size M s *°* of each symmetry sector for a 
fixed total sample size M = M s = 500. Again 

we observe, that the distribution of symmetry sectors is 
broad at high temperatures and becomes narrow for large 
values of /?. 

Figure 6 presents results for the thermal energy of 
the isotropic Heisenberg chain (3) at zero magnetic field, 
where we can exploit the non-Abelian SU(2) symmetry of 
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Figure 6. SU(2)-SYMETTS sampling of the thermal energy 
for the isotropic Heisenberg chain with h = 0, using the same 
layout as Fig. 5. 


the model. The layout of the panels and the parameters 
are chosen in accordance with those in Fig. 5 above. In¬ 
stead of S£ ot sectors, the SYMETTS are now categorized 
in terms of the total spin Stot of each SU(2) multiplet 
in the sample. Again, the upper and middle rows dis¬ 
play the average energy of each subsample correspond¬ 
ing to a fixed Stot ■ The resulting values for ( E)p deter¬ 
mined by SYMETTS [Eq. (8)] are benchmarked against 
METTS calculations [Eq. (1)] and quasi-exact purifica¬ 
tion data. We find that the overall ensemble average of 
all SYMETTS subsamples gives a good approximation of 
the thermal energy of the state also for the non-Abelian 
sampling routine. 

As for the Abelian case, the distribution of different 
multiplets shown in the last row becomes more narrow 
towards lower temperatures. Whereas the majority of 
states at /3 = 20 belong to the multiplets S to t = 0,1, 
these sectors deplete for higher temperatures and the 
maximum moves towards Stot = 2 for /3 = 4. 


longer time scales. Here SYMETTS offers a great ad¬ 
vantage over the existing METTS approach, since the 
symmetry implementation strongly increases the numeri¬ 
cal efficiency during the real-time evolution. In addition, 
the overhead cost of generating both a symmetric and 
non-symmetric sample in the SYMETTS sampling (see 
Sec. Ill C) can be ignored in almost every case, as it only 
accounts for a very small fraction [0(1CG 3 )] of the total 
computational time. The achievable efficiency gains are 
completely analogous to the exploitation of symmetries 
in other MPS applications, such as ground state DMRG, 
tDMRG or iTEBD. For example, the direct implementa¬ 
tion of the Abelian U(l) symmetries in spin models can 
already speed up calculations by about a factor of up 
to 10. 3 '-'’ 1 - 52 j] ven l ar ger benefits can be achieved when 
studying models with multiple Abelian or non-Abelian 
symmetries. 

To simulate a response function using real-time evo¬ 
lution, we follow Ref. [34] and compute for every 
SYMETTS \(j) q ) in our sample the expectation value 

m q \ e ^]B[e-^C\^ q )] (14) 

by carrying out two independent real-time evolutions, 
\ip q (t)} = e~ iHt C\4> q ) and \</> q (t)) = e~ lHt \(j) q ) using 
standard tDMRG. Equation 14 can be evaluated at any 
intermediate time-step t by calculating the overlap for 
the operator (ip q (t)\B\<p q (t)), as illustrated in Fig. 7. In 
the end, we take the sample average to obtain a result 
for the finite-temperature response function. 

In principle, there are other options for calculating the 
real-time evolution of response functions. 34 In this work, 
we restricted ourselves to the scheme outlined above, as 
it requires only two tDMRG simulations per sample state 
to access all intermediate time-steps up to the maximally 
reached timescale t max . 

Instead of studying real-time response functions, here 
we consider their Fourier transforms, i.e., spectral func¬ 
tions. More particularly, we focus on dynamic spin struc¬ 
ture factors S a P(oj,k), which are the Fourier transform 


B. Dynamic observables: Dynamic spin structure 
factor 

Whereas SYMETTS does not offer any significant com¬ 
putational advantage over the original formulation for 
computing static observables, its potential is enormous 
for the calculation of dynamic quantities, such as re¬ 
sponse functions of the form 

Af 6 {t) = (B(t)C)p , with B(t) = e iftt Be~ i6t . (13) 

For such problems, generating the ensemble states repre¬ 
sents only a negligible part of the total computational 
costs. Most computational effort has to be put into 
the real-time evolution of each state in the sample, as 
the linearly growing entanglement requires an exponen¬ 
tial increase of the bond dimension of the MPS towards 


[(cf> q \e iftt ]B[e- i6t C\^ q )] 

/ \ 

411(3 tDMRG 


. / 'A 

e~ lHtl \c/> q ) -*[...] 5 [...]*- e~ i ^ tl C\4> q ) 

I I 

tDMRG tDMRG 

+ I 

e ~ lkt ^ q ) e~ ikt *6\<l> q ) 

I i 

tDMRG tDMRG 

i I 


Figure 7. Schematic illustration of the dynamic calculation 
with SYMETTS. 34 
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of dynamical spin correlation functions. These quantities 
are of particular experimental relevance, as they can be 
directly accessed by inelastic neutron scattering experi¬ 
ment. For a benchmark, we compute the dynamic spin 
structure factor of the XXZ Heisenberg model with open 
boundaries: 


- 1 *"$7r [ * 

ij ' ' 

(i5) 

To this end, we define the spin-wave operator = 

V^lEj=i sin (^) 4 ? and evaluate via 

Eq. (14) for a number of intermediate points up to some 
maximum time t max . Then, we perform a Fourier trans¬ 
form to frequency space, including a Gaussian broad¬ 
ening exp[—4(t/t max ) 2 ] in the integral in Eq. (15) to 
remove artificial oscillations, which are caused by the 
finite cut-off of the real-time evolution . 21 This means 
that the exact spectral features are convolved with a 
Gaussian exp[—w 2 /(2W 2 )], with a frequency resolution 
W = 2\/2i“ 1 x . In some cases, linear prediction can be 
used to avoid the artificial broadening and extract more 
spectral information from the time series . 22,23 However, 
we found that linear prediction is not reliable in our study 
of the generalized diamond chain (see Sec. VB). Hence, 
we refrain from employing linear prediction in this work. 

In a first study, we employ our U(1)-SYMETTS ap¬ 
proach to extract the dynamic spin structure factor in the 
limit of A = 0. In this case, the XXZ model can be solved 
exactly by mapping the system by a Jordan-Wigner 
transformation to non-interacting spinless fermions . 53,54 
This allows us to exactly evaluate the spin correlation 
functions (S?(t)Sj (0)) for arbitrary times and obtain the 
dynamic spin structure factor by Fourier transformation 
for direct comparison to the SYMETTS data. 

Fig. 8 (a) displays the real-time evolution of S zz (7 t/4 , t) 
for (3 = 4 up to f max = 30, with the thin lines correspond¬ 
ing to individual realization of particular SYMETTS 
states and the thick red line denoting the ensemble av¬ 
erage. The sample averages are collected for different 
temperatures in Fig. 8 (b). After the real-time evolution, 
we perform a Fourier transform to obtain the dynamic 
spin structure factor as a function of frequency, as shown 
in Figs. 8 (c,d) for k « n/2 and 7 r/ 4 . We find excellent 
agreement with exact results (dashed lines). A prereq¬ 
uisite for agreement of this quality is that the statistical 
sampling error, and hence the temperature, is sufficiently 
small; for the sample size of M = 300 used in Fig. 8 the 
relative error, defined as 

\/J d uj[S(k,u}) - S exact (k,u})] 2 
< 55 = A- -, (16) 

yj f dw5 exact (£:, w ) 2 

varies between SS ~ 1 % for /3 = 4 and SS « 0.3% for 
/3 = 20. We note that the error is approximately pro¬ 
portional to the temperature, which indicates that the 



J 

\ =E 

(C) 




co Ik 


Figure 8. tDMRG calculation for the dynamic spin structure 
factor of the XX model using M = 300 ensemble states. Panel 
(a) illustrates the time evolution of individual SYMETTS 
(thin lines) at /? = 4 used to calculate S zz (n/4, t) by tak¬ 
ing the ensemble average (thick red line), (b) Displays the 
SYMETTS ensemble average for various inverse tempera¬ 
tures, which is then used to compute the dynamic spin struc¬ 
ture factor in frequency space. Panels (c) and (d) show the 
frequency data obtained from Fourier transform for k ~ 7t/4 
and n/2. For all considered inverse temperatures, we find 
excellent agreement with the exact result (dashed lines). 


dominant contribution is given by the statistical error of 
the ensemble, which scales as ~ T /\[M. 

Next, we consider an isotropic coupling A = 1, which 
allows us to compute S(k,uj ) using SU(2)-SYMETTS, 
since the XXZ Hamiltonian (3) features the full spin sym¬ 
metry in this case. Figures 9(a) and 10(b) show the 
results for k « 37 r /4 in time and frequency space, re¬ 
spectively, in comparison to purification calculations us¬ 
ing the same truncation criterion (black dashed lines). 
The maximum time f max varies for different tempera¬ 
tures, since we stopped the SYMETTS calculations when 
a threshold of 1000 states was exceeded by the bond di- 


P=50 - P=20 - p=8 - p=4 



t cd/jc 

Figure 9. SU(2)-SYMETTS calculation (solid lines) for (a) 
S(3n/4,t) and (b) S(3n/ A,ui) of the isotropic Heisenberg 
chain using M = 300 ensemble states. For all considered 
inverse temperatures, we find very good agreement with data 
obtained from matrix-product purification (dashed lines). 
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SYMETTS 


(3=50 (3=20 (3=8 (3=4 

SU(2) SYMETTS 


Purification 



(3=50 (3=20 p=8 (3=4 

U(1) SYMETTS 



Figure 10. (a) Average multiplet bond dimension D* and 

(b) the corresponding bond dimension D of SU(2)-SYMETTS 
during the calculation of ^(3^/4, t ) (solid lines) in comparison 
to the purified density matrix (dashed lines). In both cases, 
we keep all singular values > 10~ 4 during real-time evolution. 
For /3 = 4, 8 we include a backwards time-evolution of the 
auxiliary bonds of the purified density matrix, as this leads 
to a reduction of the entanglement growth. 24,55 (c) Average 
bond dimension D of U(1)-SYMETTS sample during the cal¬ 
culation of S Z2 (37 t/4 , t) for comparison. 


mension D 1 which was determined adaptively by keeping 
all singular values > 10 -4 (see Fig. 10 and upper panel 
of Table I for values of f max ). Again, we find excellent 
agreement for the considered temperature range. 

In this context, we briefly discuss the intriguing ques¬ 
tion whether SYMETTS can reach longer time scales 
than purification in certain limits. To this end, we study 
the growth of entanglement during the real-time evolu¬ 
tion, which manifests itself in the growing bond dimen¬ 
sion of both the average SYMETTS as well as the purified 
density matrix. Figures 10(a) and 10(b) present SU(2) 
data for the average maximum multiplet bond dimen¬ 
sion D* and the corresponding states space dimension 
D , respectively. We find that an average SYMETTS re¬ 
quires significantly less numerical resources at /3 = 20, 50. 
For such low temperatures, SYMETTS certainly allows 
to access longer time scales than purification when fix¬ 
ing the numerically feasible bond dimension to an upper 
cutoff. Note that due to the presence of the statisti¬ 
cal error, this does not imply that SYMETTS is gen¬ 
erally more accurate than purification when fixing the 
total computation time and judging accuracy by com¬ 
paring to quasi-exact calculations, as done in Ref. [34], 
However, if one does not insist to push the statistical er¬ 
ror towards the order of the truncation error, and more¬ 
over takes into account parallelizability, SYMETTS of¬ 
fers much potential towards the dynamical description 
of low-dimensional systems at low temperature. This is 
already illustrated by the calculations in this sections, 
demonstrating that it is possible to extract the dynamic 
structure factor with high accuracy using a sample size 
of only a few hundred states. 

On the other hand, SYMETTS is limited to small f max 
at high temperatures. Particularly at /? = 4, a single 


SU(2)-SYMETTS on average requires larger bond dimen¬ 
sions D than the purified density matrix! This can be 
attributed to the intrinsic structure of the SU(2) symme¬ 
try product states. Although their multiplet dimension 
D* is strictly unity at infinite temperature, the SPS al¬ 
ready contain some entanglement due to the presence 
of non-trivial multiplet sectors with internal structure, 
which lead to a state space dimension D > 1. Because 
of thermal fluctuations, these “excited” multiplet sectors 
appear more frequently at high temperatures, which is 
illustrated by the comparison D* and D of the corre¬ 
sponding SPS samples in the lower panel of Table I for 
different temperatures. With these non-trivial multiplets 
being present in the SPS, the subsequent imaginary- and 
real-time evolution obviously also induces more entan¬ 
glement. This explains why 5 of a SU(2)-SYMETTS 
exceeds the bond dimension necessary to represent the 
purified density matrix already before starting the real¬ 
time evolution at /3 = 4 [c.f Fig. 10(b)]. 

This issue is not present in the context of U(l)- 
SYMETTS, where the SPS does not contain any intrinsic 
entanglement at infinite temperature and thus can still 
be considered as a classical product state. Thus, the ini¬ 
tial D is strictly smaller than the bond dimension of the 
purified density matrix for all temperatures, as shown 
inFig. 10(c). Moreover, the increase of D at high and 
intermediate temperatures is slightly less severe than in 
the SU(2) calculations. 

We conclude from this analysis that it is possible for 
SYMETTS to exploit both Abelian and non-Abelian 
symmetries. In combination with tDMRG, it represents 
a valuable alternative to computing spectral functions, 
particularly for low temperatures. In addition, we find 
that Abelian SYMETTS have favorable entanglement 
properties over their non-Abelian counterparts at high 
temperatures. Hence, one should refrain from exploit¬ 
ing non-Abelian SYMETTS in these cases and switch to 
U(1)-SYMETTS or matrix-product purification. 


^max 

II 

OO 

II 

13 = 20 

/3 = 50 

SYMETTS SU (2) 

11.2 

21.4 

> 35 

> 35 

PurificatioiiRTT^i 

18 

22 

26.1 

30.6 

SYMETTS U(1) 

19.8 

31.2 

> 35 

> 35 

bond dim. 

II 

OO 

II 

/3 = 20 

O 

II 

^SU(2) 

1 

1 

1 

1 

D SU(2) 

4.11 

3.08 

2.31 

1.96 

max[L> su(2 )i 

16 

12 

7 

5 

D V(1) 

1 

1 

1 

1 


Table I. Upper panel: maximum time t max reached in the 
simulations shown in Figs. 9 and 10. Lower panel: aver¬ 
age multiplet bond dimension D* and corresponding average 
and maximum bond dimension for SU(2) and U(l) symmetry 
product states for various temperatures. 

We point out that the METTS algorithm in princi¬ 
ple can also exploit time-translational invariance in or¬ 
der to reformulate the response function in terms of 
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(B(t/2)C(—t/2)), which effectively doubles the maxi¬ 
mum reachable time scale tmax- 25 ’ 2 '’ 34 Ideally, the time- 
evolution is then carried out in the Heisenberg picture by 
evolving B and C directly in terms of matrix-product- 
operators (MPO ), 56 so that it still requires only two 
tDMRG simulations to access all intermediate time-steps. 
We note that working in the Heisenberg picture is gen¬ 
erally considered to be suboptimal for matrix product 
purification. 2. However, it seems more appealing for the 
METTS framework as one could carry out the real-time 
evolution only once for the MPO and compute the re¬ 
sponse function by calculating the overlap of the time- 
evolved MPO and the METTS sample. Thus the time- 
evolved MPO could be recycled for arbitrary tempera¬ 
tures. In general, this would imply that the maximum 
reachable time scale is set by the real-time evolution of 
the operators. In the pure-state formulation, t max would 
then be temperature independent, as temperature only 
enters through the calculation of the overlaps with the 
METTS sample. Naturally, this idea enormously profits 
from the inclusion of symmetries into the METTS lan¬ 
guage presented here, but is beyond the scope of this 
paper and will be discussed elsewhere. 

Finally, we remark that we have also explored the pos¬ 
sibility of combining SYMETTS with a Chebyshev ex¬ 
pansion to directly compute spectra in frequency space. 
However, this approach is computationally more expen¬ 
sive due to technical reasons and therefore not recom- 
mendable (see Appendix A for details). 


V. GENERALIZED DIAMOND CHAIN MODEL 
FOR AZURITE 

In the following, we demonstrate the efficiency of 
SYMETTS by studying a more complicated spin-chain 
model of direct experimental relevance. We focus on 
the natural mineral azurite Cu 3 (C 03 ) 2 ( 0 H) 2 , which 
has attracted much attention due to the discovery of 
a plateau at | in the magnetization curve at low 
temperatures . 4,5,39-49 Some authors proposed that the 
magnetic properties of this material are well described by 
a spin - 4 diamond chain formed by the copper atoms with 
purely antiferromagnetic exchange couplings . 42,45,57 Oth¬ 
ers suggested a dominant ferromagnetic coupling 4,43,44 
and the importance of interchain coupling , 46 yet none of 
them were able to derive a microscopic model for azurite 
that is able to fully characterize its complex magnetic 
properties. 

Employing a combination of first-principle methods, 
exact diagonalization, and DMRG, Ref. [39] recently 
derived a full three-dimensional model which can be 
mapped to an effective one-dimensional system, namely 
a generalized diamond chain model with purely antiferro¬ 
magnetic couplings, illustrated in Fig. 11(a). One third 
of the Cu spins (dark blue balls) forms weakly coupled 
monomers (dashed horizontal lines), whereas the other 
two thirds (light blue balls) form strongly coupled dimer 


singlets (heavy vertical lines). The dominant energy scale 
is determined by the dimer-dimer coupling J 2 . In addi¬ 
tion, there are nearest- and third-nearest neighbor dimer- 
monomer exchange J\ and J 3 as well as the monomer- 
monomer coupling J m . More precisely, the full Hamilto¬ 
nian of the generalized diamond chain is defined as 


H 0 


N/3 

m,j • (Sdl,j+1 + Sd'2,j) + J2&dl,j ' Sd2,j 

3= 1 

-\-J 3 SraJ ‘ ( &dl,j H" &d2,j-\- 1) H - 1 


N/3 


~9HbH ^ [<! 


'dl,j 


Shj + S Z mJ 


3 = 1 


(17) 


with external magnetic field H. Bohr magneton /j ,b and 
gyromagnetic ration g = 2.06. 58 N labels the total num¬ 
ber of Cu spins in the system, the number of unit cells 
is therefore given by N/3. Note that this model features 
a U(l) sp in symmetry for finite values of H, which we ex¬ 
ploit in our SYMETTS calculations. The value of the 
couplings has been determined by DFT calculations and 
small refinements using experimental data, leading to 


Ji = 15.51 K, J 2 = 33 K, 

J 3 = 6.93 K, J m = 4.62 K. (18) 

Based on this system, the authors of Ref. [39] managed to 
derive a full microscopic picture for azurite, that is able 
to explain a wide number of experimental results. Addi¬ 
tional support for the validity of this model is given by 
Ref. [40], which explores further aspects such as magne¬ 
tocaloric properties and excitation spectrum. Although 
Refs. [39] and [40] also present some selected results 
for the dynamic spin structure factor using dynamical 
DMRG (DDMRG), their resolution in the energy oj and 
the momentum transfer k is limited, since DDMRG is 
numerically expensive and requires separate calculations 
for each ui. Moreover, their results are restricted to zero 
temperature. 

We will now illustrate the power of SYMETTS by cal¬ 
culating the excitation spectra in the plateau phase and 
analyze the influence of magnetic field and temperature 
on the excitation branches. 


A. Magnetization plateau 

The most striking feature of azurite is the plateau at 
I in the magnetization curve as a function of a magnetic 
field applied along the b axis of the crystal lattice . 42 This 
property can be nicely captured by the generalized dia¬ 
mond chain model, as already demonstrated in Ref. [39] 
by employing a direct comparison of the magnetization 
obtained from ground-state DMRG calculations with ex¬ 
perimental data. Performing DMRG on an open chain 
with N = 90 spins, we obtain the magnetization plateau 
at T = 0 shown in Fig. 11(b) (black line). For bet¬ 
ter comparison, we use experimental units in the rest of 
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unit cell j 




temperatures, the quantum mechanical properties of the 
system are almost erased and the magnetization curve 
becomes a linear function of the magnetic field. We 
note that in the plateau phase, the effect of temperature 
on the magnetization depends strongly on the specific 
field strength. For values of H significantly smaller 
than 20 T the magnetization strongly decreases with 
increasing temperatures. In this case, we expect the 
change in magnetization to be predominantly caused 
by the monomers, which have to vacate their fully 
polarized state. On the other hand, the monomers 
cannot contribute to the thermal increase of the total 
magnetization for H > 20 T as they are already fully 
polarized on the plateau. Here the thermal fluctuations 
should predominantly excite the dimers by breaking up 
their singlet structure. We expect this to be reflected in 
the excitation spectra on the plateau. We study these 
next by means of the transverse dynamic spin structure 
factor, for the three values H = 14, 20, 27 T indicated 
by the vertical dashed lines in Fig. 11(b). 


Figure 11. (a) Illustration of the generalized diamond 

chain model with the antiferromagnetic exchange couplings 
Ji, Ji, J 3 , and Jm ■ One unit cell of the system is highlighted 
by the gray area, (b) Dependence of the magnetization on an 
external magnetic field H, which we calculated by employ¬ 
ing DMRG and SYMETTS at zero and finite temperature, 
respectively. In all calculations, we keep every singular value 
larger than the truncation threshold s to1 = 10 -5 and use a 
sample size of M = 1000. For a system of N = 90 spins in to¬ 
tal, the emergence of the | plateau can be observed for fields 
in the range of Hi iC ^ H Si H U:C . At finite temperatures, the 
plateau is washed out and the magnetization curve becomes 
a linear function of H. The vertical dashed lines indicate the 
parameter choices for our dynamical SYMETTS calculations 
inSec. VB. 


this section. We find that the plateau phase is bounded 
by a lower and upper critical field, H/ iC ss 9.8 T and 
H u c ~ 31.0 T. Note that the small intermediate steps 
for H < Hi c and H = 28.5 T are artifacts caused by 
finite-size effects of the chain. 

The plateau can be explained by a very intuitive 
argument.' 59 The dominant dimer-dimer exchange cou¬ 
pling J 2 forces the dimer spins into a singlet state, 
whereas the monomer spins are only weakly coupled by 
J m . Therefore, the monomer spins polarize first for a 
finite magnetic field, whereas the dimers remain in the 
singlet state for a considerable interval of H. Only at 
large fields H > H u c , the dimers are arranged in a po¬ 
larized state. Thus only | of the total spins is aligned 
in direction of the field at intermediate fields strengths 
Hi c < H < H Ujc , leading to the emergence of the | 
magnetization plateau. 

Introducing thermal fluctuations by employing 
SYMETTS, we observe that the plateau is gradually 
washed out with increasing temperature. At high 


B. Transverse dynamic spin structure factor 


We employ U(1)-SYMETTS to compute the transverse 
dynamic spin structure factor of the generalized diamond 
chain model, which can directly be measured by neutron 
scattering experiments. Following Ref. [40], the dynamic 
spin structure factor is defined as 


m,n 


^ik(Ri — Rj) 


d t e iut (S?{t)S?) 


. ( 19) 

Note that it is important to use the precise positions 
Ri of the Cu spins in azurite 59 and the experimentally 
chosen momentum direction in order to make the data 
comparable to the experiment in Ref. [4]. 

We perform all calculations for an open chain of N = 
90 spins, which allows an accurate resolution of the mo¬ 
mentum transfer k along the chain direction. For each 
k, we average over a SYMETTS sample of 300 states ex¬ 
ploiting the U(l)spi n symmetry of the model. In compar¬ 
ison to non-symmetric METTS, U(1)-SYMETTS yields 
a reduction of CPU time by a factor between 4 and 10 
for the parameters considered here (see Appendix B for a 
more detailed assessment). Using a second-order Trotter 
decomposition, we set the time step rg = Td yn = O.O 5 J. 7 1 
and truncation error s]g o1 = 10 -5 , = 5 x 10 -4 in the 

imaginary- and real-time evolution, respectively. We stop 
the real-time evolution at t max = 50 J^ 1 and checked 
that calculations are not impaired by finite-size reflec¬ 
tions on this time scale. This setup leads to a maximum 
bond dimension D < 600 at f max for all the time-evolved 
SYMETTS considered. To minimize finite-time effects, 
we here use a Gaussian broadening to perform the Fourier 
transform, leading to a frequency resolution W ~ 0.16 
meV. As an alternative route, we had also tested linear 
prediction, but found that for this model its results were 
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H= 14 T H= 20 T H= 27 T 



0 jt/5 2 jt/ 5 3jt/5 4it/5 jt 0 jt/5 2ji/ 5 3it/5 4 jt/ 5 Jt 0 jt/5 2 jt/5 3ji/5 4 jt/5 jt 


k k k 

Figure 12. Transverse dynamic spin structure factor S xx (k,u}) of a generalized diamond chain model for azurite. The intensity 
is displayed in arbitrary units. Each column indicates a different magnetic held strength, H = 14, 20, 27 T, corresponding to 
distinct points in the | plateau phase. The ground state spectra are displayed in panels (a)-(c) in the first row, whereas the 
panels (d)-(i) show finite-temperature results obtained using tDMRG in combination with U(1)-SYMETTS (for details see 
text). 


very sensitive to changes of the regularization parame¬ 
ter and the statistical window on the given time scale. 
Therefore while in principle, after significant further fine 
tuning, linear prediction may allow a systematic extrap¬ 
olation to longer time scales to enhance spectral resolu¬ 
tion, we did not further pursue this route. 

Our results are displayed in Fig. 12. Each column indi¬ 
cates a different magnetic held strength H = 14,20,27 T 
and each row corresponds to a different temperature 
T = 0,4.125,8.25 K. 

Let us first note that Fig. 12(a) at zero temperature 
and H = 14 T nicely reproduces all features of the 
DDMRG data used in Ref. [39] for a direct comparison 
with the experiment in Ref. [4]. We observe a gapped 
system with a low-energy band dispersing along fc, cor¬ 
responding to the monomer excitations, and a dimer 
branch at higher energies, whose dispersion is weakened 
by the competition of J\ and J 2 . 39 The spectral weight 
in both branches is mainly distributed around k = ir. 
Moreover, we find an additional excitation branch at 
w > 3 meV with only small spectral weight and almost 
no dispersion. In the experiment, this branch is shifted 
towards higher energies (by ~ 1 meV). 4 

Increasing the magnetic field has an effect on the po¬ 
sition of both the monomer and dimer band (but not on 
their dispersion), which can be understood easily in an 
intuitive picture. As discussed in Sec. V A, the monomers 
are fully polarized in the entire plateau phase. Hence, ex¬ 
citing a monomer spin becomes increasingly expensive for 


larger magnetic fields, because a spin-flip is penalized by 
the additional Zeeman energy. Comparing the position 
of the monomer branch in Figs. 12(a) and 12(c), the shift 
towards higher energies at H = 27 T is fully captured by 
the change in the Zeeman term « 1.6 meV. 

The magnetic field has the reversed effect on the dimer 
band, which is shifted to lower energies. Again the effect 
can be understood using the same line of arguments. Ex¬ 
citing a dimer singlet results in the break off of the singlet 
structure, allowing the dimer spins to polarize in the di¬ 
rection of H. At larger field strength, each excited dimer 
spin is therefore rewarded by a factor of (l/2)g/rsA7L 
from the Zeeman term. This fully accounts for the shift 
of the dimer branch to lower energies in Figs. 12(a) and 
12(c). At H = 20 T, the system is approximately probed 
in the middle of the plateau phase (c.f. Fig. 11). At this 
point, the band gap reaches a maximum, A E ss 1 meV, 
as the monomer branch has already moved to rather high 
energies while the dimer band is about to cross it, as il¬ 
lustrated in Fig. 12(b). 

Based on this discussion, we can confirm the very dis¬ 
tinct effects of temperature on the different points at the 
I plateau and put the arguments given in Sec. V A on 
solid ground. For regions of the plateau where the mag¬ 
netization decreases at finite temperature, thermal fluc¬ 
tuations primarily excite the monomers as this is energet¬ 
ically favorable. On the other hand, the thermal increase 
of the magnetization for larger magnetic fields observed 
in Fig. 11(b) can be understood in terms of the lowering 
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Figure 13. Transverse dynamic spin structure factor S xx (k, cj) at k as 47r/5 of a generalized diamond chain model for azurite for 
three different field strengths in the plateau phase. At the edges of the plateau phase [(a) and (c)], the large peaks indicating 
the monomer and dimer branches are already washed out at intermediate temperatures. Thermal fluctuations redistribute 
spectral weight in between the two excitation peaks in these cases. In the middle of the plateau phase (b) thermal broadening 
is almost not present at T = 4.125 K. Only at high temperatures is the height of the combined peak of monomer and dimer 
excitations significantly reduced. 


of the dimer excitation energy due to the additional re¬ 
wards in Zeemann energy, which has the opposite effect 
on the monomer band. 

Fig. 12(d-i) displays the evolution of the spin excita¬ 
tions at finite temperature. The thermal broadening ef¬ 
fects are strongly visible at H = 14 and 27 T, where 
the low energy bands are strongly smeared out even at 
intermediate temperatures. This is expected from the 
comparatively small band-gap at T = 0 and Fig. 11(b), 
which shows strong effects of temperature on the mag¬ 
netization in this regime. In contrast, thermal fluctua¬ 
tions have a much weaker effect at H = 20 T, where the 
band gap is maximal. Indeed, comparing Figs. 12(b) and 
12 (e), we see almost no difference in the distribution of 
the spectral weight. Only Fig. 12(h) shows some thermal 
broadening, yet no new features arise in the spectrum. 
Again, this is in good agreement with the robustness of 
the magnetization for finite temperature in the middle 
of the plateau, as illustrated in Fig. 12(b). These fea¬ 
tures become even more prominent when studying cross 
sections of Fig. 12, i.e., the spin excitations for a spe¬ 
cific momentum value. These are displayed in Fig. 13 
for k ~ 4/57T. Again, we observe that the large peaks 
indicating the monomer and dimer branches are already 
washed out at intermediate temperatures at the edges of 
the plateau phase [Fig. 13(a) and 13(c)]. In both cases, 
thermal fluctuations strongly redistribute spectral weight 
in-between the two excitation peaks. In contrast, the 
height of the combined monomer and dimer excitation 
peaks in the middle of the plateau phase is significantly 
reduced only at high temperatures [Fig. 13(b)]. 

To conclude, our finite-temperature study of the spec¬ 
tra of the generalized diamond chain model for azurite 
fits in nicely with previous work 39,40 and provides new 
insight in the plateau phase. We observe a crossing of 
the monomer and dimer branches with increasing mag¬ 
netic held, which can very intuitively explain the effects 
of finite temperature on the plateau phase. Testing these 
features in neutron scattering experiments would provide 


additional information on the validity of the microscopic 
model for azurite. Such a study would be particularly 
enlightening in the context of the results provided by 
Ref. [5] , which showed discrepancies of using an isotropic 
spin model to describe azurite in the regime of H < H i,c 
i.e., for Helds below the plateau phase. 


VI. CONCLUSION 

In this work, we have introduced an intuitive and easily 
implemented extension of the minimally entangled typi¬ 
cal thermal state approach of Ref. [30] , which allowed us 
to generate a METTS sample of symmetry eigenstates. 
We explicitly showed how to construct such a SYMETTS 
ensemble exploiting both the Abelian U(l) sp i n and non- 
Abelian SU(2) sp i n symmetry of spin-^ Heisenberg chains, 
without introducing strong autocorrelation effects in be¬ 
tween the ensemble states. 

Whereas SYMETTS does not improve the numerical 
efficiency when calculating static observables as com¬ 
pared to METTS, the benefits of using symmetries come 
fully to the fore when calculating more complex dynamic 
quantities such as response functions. Here, most compu¬ 
tational effort has to be put into the real-time evolution 
of each state in the sample and the gains of explicitly ex¬ 
ploiting symmetries in the MPS simulations is enormous. 
We checked the validity of our approach for the dynamic 
spin structure factors of the XX and XXX Heisenberg 
chains and found that SYMETTS in principle is able to 
reach longer time scales than purification at low temper¬ 
atures. 

Moreover, we applied SYMETTS to study the finite- 
temperature excitation spectra of a generalized dia¬ 
mond chain model for the natural mineral azurite 
Cu 3 (C 0 3 ) 2 (OH) 2 . Focusing on the plateau phase of the 
system, we found very distinct effects of temperature on 
the different points at the | plateau, which are caused 
by the Zeeman term shifting the dimer and monomer 
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branches in opposite directions. Our results fit in nicely 
with previous work 39,40 and provide new insight in the 
plateau phase. 

Interesting questions for future work involve the treat¬ 
ment of fermionic systems, where the symmetric ensem¬ 
ble states could be formulated in terms of a combina¬ 
tion of SU(2) charge and SU(2) spin symmetries or their 
Abelian counterparts. For example, SYMETTS could 
be employed to study finite-temperature density profiles 
in interacting quantum-point contacts. 60 In this context, 
it would be particularly interesting to further explore the 
possibility of combining a real-time evolution to an MPO 
with local support in the Heisenberg picture, as briefly 
described at the end of Sec. IV B. In principle, this would 
simplify combining METTS with the concept of time- 
translational invariance 2 ' to double the maximum reach¬ 
able time scale and could be a generally more efficient 
approach for finite-temperature response functions at low 
temperatures. 

Finally we note as an outlook that the SYMETTS algo¬ 
rithm may also be entirely based within symmetry eigen¬ 
states, in that the non-symmetric sampling as described 
in this paper is fully replaced by Metropolis sampling. 
Based on the weights P q above, necessarily, this must 
also include a proposal distribution to switch to neigh¬ 
boring symmetry sectors. To minimize rejection prob¬ 
ability, this random walk towards neighboring symme¬ 
try sectors can be chosen temperature dependent. By 
definition, the Metropolis sampling also guarantees de¬ 
tailed balance. And, by rejecting certain higher-energy 
states, this may lead to reduced spread and hence en¬ 
hanced convergence of computed observables. In this 
formulation, SYMETTS would also provide benefits for 
the calculation of static properties and might allow the 
finite-temperature treatment of 2D clusters. 
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Appendix A: Chebyshev expansion and METTS 

Chebyshev expansion techniques have been success¬ 
fully established as an alternative approach for the com¬ 
putation of spectral functions in the context of kernel 
polynomial methods. 61 More recently, Ref. [29] intro¬ 
duced the Chebyshev expansion in the MPS formalism 
(CheMPS) to determine spectral properties at zero tem¬ 
perature. Based on this work, CheMPS has been applied 
to determine signatures of the Majorana fermion in the 


interacting Kitaev model, 62 in the context of the inter¬ 
acting resonating level model, 63 and as impurity solver 
for single- and two-band DMFT calculations in combi¬ 
nation with linear prediction. 64,65 In addition, CheMPS 
has been expanded towards finite-temperature calcula¬ 
tions using a Liouvillian in a matrix-product purification 
framework. 28 

The question as to whether CheMPS is the most effi¬ 
cient method for computing spectral functions using MPS 
methods cannot be generally considered settled, as there 
is no one-to-one correspondence of CheMPS in its most 
efficient setup to real-time evolution. Nevertheless, the 
claim of Ref. [29] that CheMPS is significantly less ex¬ 
pensive than tDMRG to obtain the same spectral infor¬ 
mation can no longer be supported. 66 We have not con¬ 
ducted a systematic comparison of both approaches, but 
in our experience CheMPS and tDMRG require similar 
computational effort when aiming for the same spectral 
resolution and employing an equal truncation criterion at 
zero temperature. CheMPS, though, offers a significant 
advantage over tDMRG as it allows better control over 
the broadening procedure of the spectral data. 29,63 

In this context, it is worthwhile to explore the com¬ 
patibility of SYMETTS and CheMPS. To this end, we 
start with the Fourier transform of the response function 
in Eq. (13) 

Af d (uj) = J duj{5{id - H)BS(uj + u- H)C)p . (Al) 

To compute the response function in this form, we follow 
Ref. [61] and expand both 5-functions in terms of orthog¬ 
onal Chebyshev polynomials of the first kind T m (u) + ui) 
and T n (u>) before integrating over the frequency index to 
for every SYMETTS \<j> q ) in our sample. This “double” 
Chebyshev expansion involves Chebyshev moments of the 
type = {T m (H')BT n (H')C) p, where H' represents 
the Hamiltonian with a rescaled spectrum u' £ [—1,1] to 
ensure the convergence of the Chebyshev recursion. This 
is usually achieved by using a linear rescaling with the 
parameters a, b 


a 

The moments fi mn are determined by calculating a first 
set of Chebyshev vectors up to the desired expansion or¬ 
der Ache via the recursion relation 

I t m ) = 2H'\t rn _i) — |f m — 2 ), 
l*o) = \<t> q ), \ti) = H'\to), (A3) 

and keeping it in storage. Then, we iteratively obtain a 
second set of vectors | t n ) using a different starting vector 
|to) = C\<t> q ). For each | t n ), we compute the overlap 

fJ“mn — (tm | |tn) for ^ — 0, 1 ? Aqhe 1 • 

The sample average of the Chebyshev moments is then 
used to compute the finite-temperature response function 
in frequency space. 
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Figure 14. CheMETTS calculation for spin structure factor 
of the XX model with N = 50, s^ n = 5 x 10~ 4 , IVche = 300, 
77 = 0.1 and M = 300. For all considered inverse tempera¬ 
tures, we find excellent agreement with the exact result. How¬ 
ever, note that the required numerical resources clearly exceed 
those used in Fig. 8 , where tDMRG was employed for a system 
with twice as many spins! 


To this end, we work with a finite broadening rj instead 
of the usual kernel approach for 6(u)'+ui' — H') to remove 
the artificial “Gibbs” oscillations caused by finite expan¬ 
sion order from the spectral data . 63 Note that the broad¬ 
ening has to be performed for the 5-function contain¬ 
ing the external frequency index only as u/ is integrated 
out to obtain the final result. The finite-temperature re¬ 
sponse function in Eq. (Al) then takes the form 

Ag^{u) = - Y' /w(2 - 5 m0 ) 

H a 

m,n —0 

[ d Q'—^=Tm(u')a n (z), (A4) 

J 7TV 1 — ^ 

with z = (<*/ + u>') + irj/a and a n given by 63 

rv M = _ 2/(1 + 5 n0 ) _ 

{z) n + 1 { 1 + ^V^l/z 2 ) n ^/l - l/z 2 ' 

(A5) 

In principle, this approach represents an alternative 
to the combination of tDMRG plus Fourier transform, 
which we have applied in the main part of this work. We 
illustrate this in Fig. 14, where we used U(1)-SYMETTS 
and a double Cliebyshev expansion to compute the dy¬ 
namic spin structure factor of the XX model showing ex¬ 
cellent agreement with exact calculations (dashed lines). 
However, the Chebyshev approach in the METTS for¬ 
malism involves significantly higher computational costs 
than the real-time evolution, since, in contrast to T = 0 
CheMPS, (i) the full set of Chebyhsev vectors \t m ) has to 
be stored throughout the entire calculation, and (ii) the 
number of moments increases from JVche to N^ he , also 
squaring the number of MPS overlaps to be calculated. 

Therefore, we conclude that the combination of 
CheMPS and METTS is not a competitive alternative 
to real-time evolution for calculating spectra at finite 
temperature, since the advantage of the more controlled 
broadening procedure does not outweigh the drastically 


enhanced numerical costs involved in the double Cheby¬ 
shev expansion. 


Appendix B: Numerical efficiency of SYMETTS for 
azurite 

Here we assess the numerical performance of 
SYMETTS and the existing METTS approach on an ex¬ 
plicit example. We focus on the average cumulative CPU 
time tcpu required to carry out the real-time evolution 
of one ensemble state up to t < f max = when de¬ 

termining the dynamic spin structure factor S xx (k,co) 
of azurite in Eq. (19). For simplicity, we choose the 
same model parameters as in Fig. 13(a), namely N = 90, 
H = 14 T, and k = 4/57r, and in Figs. 15(a,b) display 
the resulting average cumulative CPU times for the two 
temperatures T = 4.125,8.25 K. Using the tDMRG setup 
described in Sec. VB, each calculation was performed on 
a single core Xeon E5-2670v2 (2.50 GHz) machine with 
4GB memory. 

The explicit implementation of the U(l) symmetry in 
the SYMETTS ensemble states clearly enhances the nu¬ 
merical efficiency, resulting in an average reduction of 
CPU time by a factor of 4 for T = 4.125 K in compari¬ 
son to a non-symmetric METTS sample [c.f. Fig. 15(a)]. 
The efficiency gain increases to a factor of almost 10 for 
T = 8.25 K in Fig. 15(b), since states with a larger bond 


T = 4.125 K 



T = 8.25 K 



10 1 - 10 * 

D 


Figure 15. (a,b) display the numerical performance of 

SYMETTS and METTS in terms of the average cumulative 
CPU time tcpu as a function of U 2 , when carrying out the 
real-time evolution to determine the dynamic spin structure 
factor S xx (k,uj) of azurite in Eq. (19) using the parameters 
N = 90, H = 14 T, and k = 4/5n. (c,d) show the average 
computation time of a single time step r as a function of D. 
We note that the ensemble states have a D « 20 at both 
temperatures after the initial imaginary-time evolution. The 
dashed lines in Fig. 15(d) are guide to the eye illustrating the 
tcpu ~ D d scaling of the CPU time for larger bond dimen¬ 
sions. The employed tDMRG setup is described in Sec. VB. 
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dimension profit even more from the exploitation of the 
spin symmetry. This is illustrated in Fig. 15(c,d), where 
we show the average computation time of a single time 
step t as a function of the average maximum bond dimen¬ 
sion D. The additional imaginary-time evolution nec¬ 
essary for the generation of the SYMETTS requires on 
average only 25 and 12 seconds for T = 4.125,8.25 K, 
respectively. Thus, the overhead costs of the generating 
the SYMETTS sample are clearly negligible compared to 
the total computation time. 

In addition to benefits in terms of memory require¬ 
ment, SYMETTS enables us to reduce the CPU time 
necessary to compute the dynamic spin structure fac¬ 


tor for various momenta and magnetic fields presented 
in Fig. 12 from 0(1O 6 ) to 0( 10 5 ) hours. Consequently, 
when running 400 CPUs in parallel, SYMETTS generates 
this data in roughly one week, whereas the same calcu¬ 
lation would require almost three months in the original 
METTS formulation. Note that the factor of 10 gained 
in numerical efficiency by implementing the U(l) spin 
symmetry in simple spin-chain models has also been re¬ 
ported in Ref. [37] and [51] in the context of iTEBD 
and tDMRG, respectively. Even larger benefits can be 
achieved when studying models with multiple Abelian or 
non-Abelian symmetries. 
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